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ABSTRACT 

We use the results of a high resolution N-body simulation to investigate the role of 
the environment on the formation and evolution of galaxy-sized halos. Starting from 
a set of constrained initial conditions, we have produced a final configuration hosting 
a double cluster in one octant and a large void extending over two octants of the sim- 
ulation box. In this paper we concentrate on gravitationally bound galaxy-sized halos 
extracted from these two regions and from a third region hosting a single, relaxed 
cluster without substructure. Exploiting the high mass resolution of our simulation 
ijTibody = 2.1 X lO^/i^^M©), wc coustruct halo samples probing more than 2 decades in 
mass, starting from a rather small mass threshold: 5 x 10^^ /^"^M© < M. We present 
results for two statistics: the relationship between 1-D velocity dispersion cTy and mass 
Mo and the probability distribution of the spin parameter P(A), and for three dif- 
ferent group finders. The a^ — Mq relationship is well reproduced by the Truncated 
Isothermal Sphere (TIS) model introduced by Shapiro et al. (1999), although the slope 
is different from the original prediction. A series of <j^ — Mq relationships for different 
values of the anisotropy parameter /3, obtained using the theoretical predictions by 
Lokas & Mamon (2001) for Navarro et al. (1996, 1997) density profiles are found to 
be only marginally consistent with the data. Using some properties of the equilibrium 
TIS models, we construct subsamples of fiducial equilibrium TIS halos from each of 
the three subregions, and we study their properties. For these halos, we do find an 
environmental dependence of their properties, in particular of the spin parameter dis- 
tribution P(A). We study in more detail the TIS model, and we find new relationships 
between the truncation radius and other structural parameters. No gravitationally 
bound halo is found having a radius larger than the critical value for gravithermal 
instability for TIS halos {vt > 34.2ro, where tq is the core radius of the TIS solution). 
We do however find a dependence of this relationship on the environment, like for the 
P(A) statistics. These facts hint at a possible role of tidal fields at determining the 
statistical properties of halos. 
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1 INTRODUCTION the Cold Dark Matter (CDM) cosmological model is repre- 
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ing" events (White & Rees 1978), whereby the subunits are 
accreted in a discontinous way, or (more hkely) by a com- 
bination of the two. In the first case the typical halo pro- 
files evolve adiabatically, while in the merging scenario each 
merging "event" will induce some transients in the char- 
acteristic properties which in turn will induce some evolu- 
tion in the typical profiles, after the subunits have been ac- 
creted and destroyed. In either case, one expects that relax- 
ation processes should drive the evolution towards a quasi- 
equilibrium state on a dynamical timescale, a state charac- 
terized by relationships among global quantities related to 
the halo, like its mass, density, velocity dispersion Gv^ and 
possibly others. Recently a considerable attention has been 
devoted to the study of one of these relationships: the radial 
dependence of the (spherically averaged) density, also known 
as the density profile. Unfortunately the density profile is a 
very difficult tool to use when trying to characterize the sta- 
tistical properties of halo populations, because the predic- 
tions of different models of halo formation differ only in the 
behaviour in the central parts, where the statistics is typi- 
cally poor. Less attention has been paid to another global 
quantity, namely the velocity dispersion and to its relation- 
ship with other global quantities, like the mass. The velocity 
dispersion enters the second-order Jeans equation, while the 
density profile is described by the zeroth-order Jeans equa- 
tion (Binney & Tremaine 1987). For this reason it contains 
different physical information than the density profile. Re- 
cently Bryan & Norman (1998) have looked at the ctv — Mq 
relationship for clusters, and they find a good agreement 
with the standard, singular isothermal sphere model as far 
as the slope of the relationship is concerned. Also Knebe 
& Miiller (1999) looked at this relationship, using a differ- 
ent code. Halo equilibrum models make predictions about 
the Gv — Mo relationship, but these are difficult to compare 
with observations, because some of the involved quantities 
(e.g. the velocity dispersion itself) are not directly deducible 
from observations. They can however be studied with N- 
body simulations, and one of the purposes of this paper will 
be to show that the av — Mq relationship can be used to 
discriminate among different halo equilibrium models. 
A second problem we will study concerns the dependence 
of halo properties on the environment within which they 
form. In both the hierarchical clustering scenarios mentioned 
above one could imagine that the properties of the halos do 
depend on the environment. For instance, the dynamics of 
the infall process could be affected by the average overden- 
sity of the environment within which the halo grows (Gunn 
1977), or by its shear (Buchert et al. 1999; Takada & Futa- 
mase 1999). Also typical quantities related to the merging, 
like the frequency of merging events, could intuitively be af- 
fected by the average density of the environment, at least for 
galaxy-sized halos forming within clusters. High resolution 
N-body simulations (e.g. Moore et al. 1999) show that most 
of the galaxies lying in the central (i.e. virialised) parts of 
clusters do not easily reach a velocity larger than the escape 
velocity: so, they are bound to the cluster for the largest 
part of their evolutionary history, and consequently form in 
an overdense environment. It would then be interesting to 
try to understand whether there are systematic differences 
between halos forming in clusters and in voids. 



Some of these issues have been recently discussed in the lit- 
erature. Lemson & Kauffmann (1999) have analysed the de- 
pendence of various statistical properties, including the spin 
probability distribution, on the environment, and found no 
evidence for any dependence apart for the extent of the mass 
spectrum. They divide their halos in groups according to the 
overdensity of the environment within which they are found, 
and show that the scatter diagrams between different quan- 
tities are indistinguishable among the different groups. In 
the present study, we follow a different strategy. We study 
a simulation obtained from a constrained set of initial con- 
ditions, in order to get a few clusters (and, in particular, 
a double cluster) and a large void within the same simu- 
lation box. We then extract our halos from three spatially 
disjoint regions: one containing a double cluster, a second 
one contaning a single cluster a third one containing the 
void. This is in some sense complementary to the procedure 
which Lemson & Kauffmann seem to have followed, because 
our halos are grouped according to the spatial distribution, 
rather than according to the overdensity, so they are grouped 
according to the environment within which they form. 
Very recently Gardner (2000) presented a study of the spin 
probability distribution for 6 different cosmological models 
and environments. He finds a difference between the dis- 
tributions of halos resulting from recent mergers and halos 
which did not experience mergers, almost independent of the 
environment within which they form. This could have signif- 
icant consequences for the construction of merger histories, 
and, ultimately, for the semi- analytical modelling of galax- 
ies. Similar results have been recently obtained by Vitvitska 
et al. (2001). 

The plan of the paper is as follows. In section 2 we describe 
the numerical setup of the simulations and the algorithm 
adopted to identify halos. In section 3 we describe the halo 
equilibrium models with which we compare the results of our 
simulation, and in section 4 we show the results of this com- 
parison and discuss their physical interpretation. Finally, in 
the conclusions we summarize our results and suggest some 
directions for future studies. 

In the following we will always assume a Q — 1 Stan- 
dard Cold Dark Matter model, with a Hubble constant 
Ho = lOO/iKms"^ Mpc~\ and h = 0.5. All lengths, unless 
explicitly stated, are assumed to be comoving. 



2 SIMULATIONS 

The simulation from which the data have been extracted 
has been performed using FLY (Becciani et al. 1997, 1998) 
a parallel, collisionless treecode optimised for Shared Mem- 
ory and/or clustered computing systems. FLY deals with 
periodic boundary conditions using a standard Ewald sum- 
mation technique (Hernquist et al. 1991). The algorithm 
adopted is the octal-tree algorithm of Barnes & Hut (Barnes 
& Hut 1986), with some modifications ("grouping" of cells 
belonging to the lists of nearby particles, Barnes 1987) dur- 
ing the phase of tree walking. These changes have a neg- 
ligible impact on the overall numerical accuracy, as shown 
elsewhere (Becciani et al. 2000), but they have a strong pos- 
itive impact on parallel peformance and scaling. 
We have performed two simulations starting from the same 
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Figure 1. Snapshot of the simulation box at the end of the sim- 
ulation. The scale of grey corresponds to density in logarithmic 
units. The large void is clearly seen in projection extending over 
the lower left octant. 



initial conditions. In both cases the underlying cosmologi- 
cal model is a Standard Cold Dark Matter (SCDM), with 
r^o = 1, cTs = 0.9. The main reason for this choice lies in the 
fact that the specific prediction for the Gv — Mq statistics we 
consider in the next sections were done for this particular 
cosmological model. We plan to extend our work to other 
cosmological models in future work. 

Each simulation used 256^ particles, and the box size was 
50/i~^Mpc, so that the mass of each particle is rupart = 
2.07 X 10^ /i~^M0. The simulations were designed to study 
the evolution of a Coma-like cluster, and for this purpose 
constrained initial conditions were prepared, changing only 
the softening length, which was fixed to e = 10 h~\pc and 
5 /i~^kpc, respectively. As far as the results presented in this 
paper are concerned, there are no differences among these 
two simulations, so for the rest of this paper we will con- 
centrate only on the simulation with the largest softening 
length, which we designate in the following as 16ML_1. 

Constrained initial conditions were prepared using the 
implementation of the constrained random field algorithm of 
Hoffman & Ribak (1991) by van de Weygaert & Bertschinger 
(1996). We took the same initial conditions adopted in one 
of the simulations from the catalogue of van Kampen & 
Katgert (1997). More specifically, we constrain the initial 



Figure 2. The region of the DOUBLE cluster. 

Table 1. Properties of the analyzed regions. 

Region x y z L Number of halos 



DOUBLE 


15 


12.5 


7.5 


10 


827 


SINGLE 


12 


-18 


4 


10 


786 


VOID 


-10 


-10. 


-10. 


20 


609 



All lengths are in h~^Mpc. From left to right, columns are as 
follows: label of the region, x,y,and z coordinate of its center, size 
of the region, total number of halos found. 

conditions to have a peak at the centre of the simulation 
box, with height 3.04cr, Gaussian smoothed at a scale of 
6h~^ Mpc, and a —2a void centered at (15,0,-10). The 
final configuration is shown in Figure 1. In order to study 
the environmental dependence of the properties of galaxy- 
sized halo populations, we selected three regions within the 
computational box, which we call DOUBLE, SINGLE and 
VOID. All the three are cubical with centers and sizes as 
specified in Table 1. The DOUBLE region hosts a double 
cluster, with two large parts in the act of merging by the 
end of the simulation (see Fig. 2). The SINGLE region hosts 
a more relaxed cluster without any apparent substructure. 
Finally, we have included in the analysis a significantly un- 
derdense region, the VOID, which is more extended than 
the former two, so as to contain enough halos to allow a 
reasonable statistics. 



2.1 Finding halos 

Various methods have been devised to extract halos from the 
outputs of N-body simulations. Some of these methods make 



© 2000 RAS, MNRAS 000, 1-13 



4 V. Antonuccio-Delogu et al. 



use only of particle positions, like the standard Friends- 
Of-Friends (herefater FOF) and the various versions of the 
Adaptive FOF, while others take into account also particle 
velocities (e.g. SKID, Governato et al. 1997) and/or envi- 
ronmental properties like local densities (like HOP, Eisen- 
stein & Hut 1998). Most of the results we will present later 
have been obtained using SKID because it selects grav- 
itationally bound groups of particles. In short, SKID first 
builds catalogues of groups using a standard FOF algorithm, 
selecting only particles lying within regions whose density is 
larger than a critical value 6crit- It then computes the escape 
velocity of each particle and discards those particles having a 
rms velocity larger than the escape velocity. This "pruning" 
procedure should then leave only those particles which are 
actually bound to the group, discarding those "background" 
particles which find themselves by accident at a given time 
within it. The initial linking length of the FOF phase deter- 
mines the approximate size of the groups we are considering. 
We assumed a linking length of 100/i~^kpc, corresponding 
to the typical size of a galaxy-sized object at the present 
epoch. The softening length for the calculation of the grav- 
itational potential was assumed to be the same as in the 
simulation, and the critical density 6crit was set equal to 
178/(1 + z), the value for nonlinear collapse in the Gunn & 
Gott collapse model, so that only particles from genuinely 
nonlinear ly collapsed shells should be included. 
Following a suggestion of the anonymous referee, we have 
also adopted two more halo finders to check the robustness 
of the results: an Adaptive FOF method devised by van 
Kampen & Katgert (1997) and a modified version of SKID 
which should avoid the problems posed by the original ver- 
sion. This particular AFOF halo finder selects only those 
haloes that are virialized, by specifically testing for virial- 
ization. Concerning the second method, we have modified 
SKID only in that part which engenders the input group list 
which is subsequently "pruned" of the non-gravitationally 
bound particles: in place of a standard FOF (as in the orig- 
inal version of SKID) we have adopted HOP(Eisenstein & 
Hut 1998) as input group finder. As we will see later, some 
of the relationships we find do depend on the group finder 
adopted, but those relationships holding for the equilibrium 
TIS halos are not affected by this. All the results we present 
in the following are for the final redshift of the simulation, 
z = 0.0047, unless otherwise stated. 



3 HALO EQUILIBRIUM MODELS 

The internal properties of halos formed by gravitational col- 
lapse can be described by looking at correlations among 
different physical quantities. The density profile has often 
been used to study the properties of relaxed, virialised halos, 
particularly since the findings by Navarro, Frenk, & White 
(1996) that this profile has a universal character when ex- 
pressed in dimensionless units. However, the density profile 
can be reliably determined only for halos having enough par- 
ticles in each shell to minimise the statistical fluctuations. 
For instance, Navarro et al. (1997) considered only 8 halos 
extracted from a low resolution simulation and re-simulated 
with a higher mass resolution. 
Typical N-body simulations on cluster scales tend to pro- 



duce a large amount of halos, whose density profile can not 
be reliably determined, because each of them contains on av- 
erage less than 10^ particles. For this reason we have chosen 
to study relationships involving global halo properties. This 
choice is not free from potential problems: systematic biases 
can be introduced by the particular group finder adopted. 
Consider for instance SKID, which works by stripping out 
gravitationally unbound particles from halos builded using 
FOF: the group catalogues so produced tend to be more bi- 
ased towards less massive halos than catalogues produced 
using FOF. We have then decided to adopt three different 
group finders, in order to be able to understand the role of 
these systematic factors. We have also considered two dif- 
ferent statistics to characterize the properties of dark mat- 
ter halos, and particularly their equilibrium properties at 
the end of the simulation: the internal 1-dimensional veloc- 
ity dispersion - mass relationship (av — Mq ) and the spin 
probability distribution P(A). Theoretical predictions con- 
cerning both of them are available in the literature. In par- 
ticular, we will compare the results from our N-body simu- 
lations with four models: the Standard Uniform Isothermal 
Sphere (SUS) (see e.g. Padmanabhan 1993, chap. 8 for a 
detailed treatment), the Truncated Isothermal Sphere mod- 
els (TIS) recently introduced by Shapiro et al. (1999), the 
"peak-patch" (PP) Montecarlo models by Bond & Myers 
(1996a), and some models derived from the Navarro et al. 
(1996, 1997) (hereafter NFW) density profiles. The first two 
models predict a cr^ — Mo relationship given by: 



^/rl/3 /-, , xl/2 7 1/3 T^ -1 

cr^=c/Mi2 (1 + ^con) ri ' Kms 



(1) 



where the subscript / = SUS or TIS, while for the PP 
model the relation is given in Bond & Myers (1996b): 



iv/rO.29 /-, I \1/2 7 1/3t^ -1 

cr^ = cppMi2 {l-\-Zcoii) h Kms 



(2) 



In the equations above M12 is the mass in units of 10^^ M© 
and Zcoii is the collapse redshift. The coefficients for these 
cases are given by 



csus,Tis,pp = (71.286, 104.69, 117.60) 



(3) 



respectively. We restrict our attention to these four models 
because the physical ingredients which enter in their formu- 
lation are very different, and encompass a sufficiently wide 
range among all the possible nonlinear collapse and virial- 
ization mechanisms. This wide choice reflects our generally 
poor level of understanding of the nonlinear physics of gravi- 
tational collapse, of its dependence on the local environment 
and on other properties like the merging history of the sub- 
structures. 

The SUS model is based on the spherical nonlinear collapse 
model (Gunn & Gott 1972). In this model the collapse to- 
wards a singularity of a spherically simmetric shell of mat- 
ter in a cosmological background is halted when its radius 
reaches a value of half the maximum expansion radius. The 
velocity dispersion is then fixed by imposing the condition of 
energy conservation, which must hold in the case of collision- 
less dark matter as that envisaged here. The TIS model also 
considers the highly idealized case of a spherically symmet- 
ric configuration, but assumes that the final, relaxed system 
is described by an isothermal, isotropic distribution func- 
tion and that the density profile is truncated at a finite 
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radius. Shapiro et al. (1999) have shown that this config- 
uration could arise from a top- hat collapse of an isolated 
spherical density perturbation if, as shown by Bertschinger 
(1985), the dimensionless region of shell crossing almost co- 
incides with the region bounded by the outer shock in a ideal 
gas accretion collapse with the same mass (in a Q = 1 CDM 
model). The truncation radius is then assumed to coincide 
with the region of shell crossing, and this allows them to 
specify completely the model. 

The peak-patch models introduced by Bond & Myers 
(1996a) are more general than the SUS model, in that they 
include a more realistic collapse model where deviations 
from spherical symmetry are taken into account. The density 
perturbation is approximated as an axisymmetric homoge- 
neous spheroid. Coupling between the deformation tensor 
and the external and internal torques are consistently taken 
into account up to a few first orders, and Montecarlo real- 
izations are used to build up catalogues of halos. These have 
been compared with the N-body simulations of Couchman 
(1991) in order to properly normalize the statistics. Note 
that eq. 2 is a best-fit relationship and holds for a range of 
halo mass (2.5 x 10^"^ < M^ < 5 x lO^^M©) much larger 
than that considered here. Nonetheless, we include it into 
our comparison because the physical model it is based on is 
significantly different from the other models we consider. 
Finally, we have considered models for the Gv — Mq statistics 
consistent with the NFW density profile, which were recently 
derived by Lokas & Mamon (2001) solving the second order 
Jeans equation: 



--7-{par)+2f3— = --— 
p dr r dr 



(4) 



where [3=1 — aQ{r)/a'^{r) quantifies the anisotropy of the 
velocity dispersion. For a NFW density profile we have: 



P°cC^9ic) 



3 s{l+csy 

In (1 + cs) 



(5) 



(6) 



pNFw{r) = pNFw{s) - 

GMviR I X 
rviR <5 

where we have defined: vvir^Mvir as the virialisation ra- 
dius and mass, respectively, s = r/rviR, c = c{M,z) is the 
concentration parameter and: g{c) = l/(ln(l + c)— c/(l + c)). 
Eq. 4 can be solved by quadrature, and the solution finite 
in the limit r ^ oo is: 



Vv 



-S,f3: 



St) 



g{c){l + cs) 



2 1-2/3 



,,2/3-3 



ln(l + cs) 



^2/3-2 



(1 + cs) 



(1+CS)3 



ds. 



(7) 



Note that we have always chosen as critical treshold for our 
group finders the virialisation overdensity {6 > 6crit = 178), 
so the quantities rviR, Mvir are the actual radius and mass 
found by the group finders for each halo. 
In order to make use of eq. 7 we have yet to specify the 
dependence of the concentration parameter on the mass at 
the final redshift: c = c{M,z = 0). We adopt the relation- 
ship provided by Bullock et al. (2001), by running their code 
CVIR for the relevant cosmological model. In the mass range 
we are interested to (5 x 10^° < M/Mq < 5 x 10^^) we find 
a power law fit: c{M,z = 0) = 472.063 x M"°-^2'^^°°^ Fi- 



(8) 



nally, in order to make a proper comparison with the quan- 
tity computed by the group finder, we evaluate the mass 
averaged velocity dispersion: 

2 _ ^^ Jo ^l,NFw{s)p{s)s'^ds 

'^"~ M(T) 

where we have defined: 

^v,NFW = CTr + C^0 = (2 - 0)crr 

and: 

M(l) = 47r / p{s)s^ds 



4 STATISTICAL PROPERTIES 



4.1 The av 



Mo relation 



In Figure 3 we show plots of the final av — Mo relationship 
for galaxy-sized halos in the DOUBLE, SINGLE and VOID 
regions, respectively, obtained using SKID. The most strik- 
ing difference is probably the different character of halos in 
the VOID region when compared with the clustered regions. 
Halos in VOID have a much smaller dispersion around the 
mean than halos in the clustered regions, and the distri- 
bution is almost symmetric with respect to the best fitting 
approximation. Halos in the DOUBLE region have a larger 
dispersion and they do show an asymmetry in the distribu- 
tion around the best fit solution, i.e. an excess of low mass 
halos at low ay. The latter point is useful to understand 
the potential systematic effects introduced by a particular 
group finder. In Figure 4 we show the ay — Mo relationship 
for the DOUBLE region obtained by using the two other 
group finders mentioned above. As is evident, the excess 
of low-mass halos is only an artifact introduced by SKID, 
which makes use of a FOF algorithm to build up an input 
list of groups. The main results of the next sections, how- 
ever, do not depend on the particular group finder adopted, 
because we will select subsamples of halos which can be re- 
garded as equilibrium TIS halos, and for them the slope of 
the av — Mo relationship is independent of the particular 
group finder initially adopted. 

A more important difference is evident from a comparison 
between clustered and void regions. Halos in VOID have a 
larger mass extent (a property already noted by Lemson & 
Kauffmann 1999) and also the slope of the av — Mo relation 
seems to be larger than for the other two cases. 
In Figure 5 we show a comparison with some theoretical pre- 
dictions for NFW density profiles. In order to apply eqs. 6, 7 
we have still to specify a relationship between the virial ra- 
dius and the corresponding mass, which enters into eq. 6. 
We do this by fitting a power law relationship to the data 
obtained from the simulations (Figure 7), which shows that 
the slope depends signficantly on the environment. As is 
evident from Figure 5, none of the models fits adequately 
the data. Note that there is only a slight difference bewteen 
isotropic (/3 = 0) and anisotropic (/3 = 0.5) models. Only if 
we allow for an unrealistically low value of the slope of the 
rviR — Mvir relationship we get a limited agreement for 
the DOUBLE cluster, but not for the other two regions. 
The slopes of the av — Mo relationship for different regions 
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Table 2. Least-square best fit parameters for cr^ relationship. 
Region a Ao; cq Method 



DOUBLE 


0.38 


0.08 


74.13 


SKID 




0.39 


0.05 


87.74 


SKID with HOPinput 




0.42 


0.07 


72.48 


AFOF 




0.42 


0.03 


82.16 


TIS selected halos 


SINGLE 


0.35 


0.04 


80.28 


SKID 




0.37 


0.04 


82.15 


SKID with HOPinput 




0.42 


0.07 


71.22 


AFOF 


VOID 


0.39 


0.04 


86.81 


SKID 




0.40 


0.04 


89.43 


SKID with HOPinput 




0.45 


0.05 


63.12 


AFOF 




0.38 


0.02 


88.60 


TIS selected halos 



a, Co are the fitting parameters of a power law fit of the form: 
(7v = coMp2 5 ^^ is the rms error associated with a. 

and using different group finders seem to be consistent with 
each other, within the errors (Table 2). None of the theo- 
retical models we are considering, however, seems to offer 
a good fit for all the cases. The TIS model seems to give a 
good fit for the DOUBLE region and for the SKID group 
finder, but when we use the modifed SKID group finder, 
which produces a sample over a more extended mass inter- 
val, we see that the original TIS model does not offer a good 
fit (Figure 4). 

Note that the rms uncertainty of ay in Figures 3-7 is 
less than 30 Km/sec, a value much lower than that found by 
Knebe & Miiller (1999) in their simulations (see their Fig. 
3) . We believe that this is a consequence of the larger mass 
and force resolution of our simulation, and also of the use of 
a larger dynamic range than adopted by previous authors. 

Generally speaking, a power-law seems to offer a good 
fit for all the three regions (although with different ranges for 
the three regions), but in order to determine the slope one 
must probably go a step further in modelling the physical 
state of these halos. In the next section we will explore in 
more detail the properties of TIS halos and we will focus our 
attention on their statistical properties. 



4.2 Comparison Avith the Truncated Isothermal 
Sphere model 

We will now consider the possibility of obtaining a rea- 
sonable fit of the av — Mo relationship by modifying the 
minimum-energy TIS model. We will then present here some 
more features of this model. 

Following Shapiro et al. (1999), the TIS solution is obtained 
by imposing a finite truncation radius rt on an isother- 
mal, spherically symmetric collisionless equilibrium config- 
uration. Shapiro et al. define a typical radius: 



ro = 



i^TTGpo) 



1/2 



(9) 



where po is the central density (TIS models are non- 
singular). Combining the Poisson and the Jeans' equilibrium 
equations, and making the hypothesis of isothermality for 
the distribution function, they obtain an equation for the 
dimensionless density (see Shapiro et al. 1999, eq. 29): 
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Figure 3. 1-D velocity dispersion versus mass for halos extracted 
from the three regions. The three fitting curves corrsond to the 
three cases considered in the text: Truncated Isothermal Sphere 
(continuous line), Bond & Myers 1996 (dotted line), Standard 
Uniform Sphere (dashed line). Note the larger mass extent for ha- 
los in the VOID. In all the cases, the slope is larger than predicted 
by models, although within clusters the statistical uncertainty is 
large. 
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Figure 4. 1-D velocity dispersion versus mass for the DOUBLE 
cluster region, halos selected using the AFOF and modified SKID 
with HOPinput group finders. The same parameters adopted for 
SKID in Figure 3 are adopted. The best fit power law for the 
plot in the upper figure has a slope: a = 0.039 ± 0.05. 



__ .^ , , _^. 



where we have used the definitions: 



(10) 



P . r 

P= — , C = — , 
po ro 

Shapiro et al. have shown that nonsingular solutions of eq. 10 
form a one-parameter family depending only on (t = rt/ro. 
The total mass is then given by: 



Mo = M(rt) 



f 

Jo 



Anp(r)r^dr — 4:7TporlM{Ct 



(11) 



where we have defined a dimensionless total mass: 



M(CO 



ro 



dr-^ 
po V^o 



10 



o 
a; 



^ li 



10 




10 



11 



ni2 



10 

group ^ 0-^ 



Figure 5. The c^^, — Mq relationship for three models based on 
the NFW density profile. Continous line: [3 = 0, slope of the 
rviR — MviR given by o; = cxmean- Dotted hue: f3 = 0.5, o; = 
Oimean- Dashcd line: f3 = 0,a = OLmean — Aq; (i.e. the la limit). 
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Figure 6. The truncation parameter Qt as a function of group's 
mass. 

We follow further Shapiro et al. (1999, eq. 41) and write 
the virial theorem for a collisionless truncated isothermal 
sphere: 



where: K = Mo{v^)/2 



(12) 



(3/2) Moa^ and Sp is the surface 
pressure term. In the Appendix we show that, after some 
simple algebra, starting from eq. 12 one obtains the following 
relation: 



GMo _ SMiCt) - piCt) i-. 



M(Ct) = $(Ct) 



(13) 



The function ^(Ct) is specified in the Appendix. 
We have already seen that a power- law fit describes well the 
ay — Mq relationship. Then from eq. 13 we deduce that also 
the truncation radius rt has a power -law dependence on 
the total mass: rt oc c{C,t)MQ~ ^ . Using the values of a from 
Table 2 we see that the slope of this relationship should lie 
within the range 0.22-0.3, i.e. it should be quite small. In 
Figure 7 we plot rt as a function of Mq. The best-fit values 
we find for the slope are inconsistent with the predictions 
from the ay — Mq relationship. 

The right hand side of equation 13 depends only on the di- 
mensionless truncation radius ^t? which in the minimum- 
energy TIS solution of Shapiro et al. (1999) should be 
fixed to Ct = 29.4. The function $(Ct) has a singularity at 
C ~ 0.97, where the denominator goes to zero (Figure 8). 
However, we are interested in the region ("t > 1, where 
the truncation radius is at least comparable to the core ra- 
dius. As is clear from Fig. 10, the function has minimum at 
C ~ 59.5. 

We will then look for solutions in the interval 2.91 < 
Ct < 59.5. Within these limits the function $(C) is mono- 



Figure 7. Relationship betwen truncation radius and mass. Data 
are fitted using a power-law relationship: rt = cM'^ . Least- 
squares fit values are: csingle = 1.0789 x ^^~^ ^olsingle = 
0.2637 ± 0.1335; cvoiD = 2.7989 x 10'^, ayoiD = 0.3816 ± 
0.0789; CDOUBLE = 7.3451 x 10'^, aoouBLE = 0.2791±0.1193. 



tonically decreasing and the solution of eq. 13 is certainly 
unique. Note that the TIS solution is unstable for (t > 34.2 
(Shapiro et al. 1999; Antonov 1962; Lynden-Befi & Wood 
1968), so our choice of the upper limit will allow us to ver- 
ify a posteriori this prediction. Equation 13 can be inverted 
w.r.t. (t given the left-hand side, because for each group 
all the quantities in the left-hand side are computed by the 
group finder. As for the truncation radius rt we adopt the 
groups' radius as computed by SKID, which coincides with 
the virial radius rviR- 

An interesting feature of the TIS solution, which is evident 
from Figure 8, is that in order to invert eq. 13 to find Ct , 
the value of the left-hand side must lie within a rather small 
range of values: 2.75 < GMo/rto-y < 4.56. As we can see 
from Figure 9 for the case of the DOUBLE region, this is 
not the case for all the halos, even for those halos closely 
verifying the ay — Mq relationship. For instance, out of 827 
halos identified by SKID within the DOUBLE region, only 
382 lie within the region for which eq. 13 can be inverted. 
For those halos for which eq. 13 can be inverted, we are able 
to determine Ct and to compare it with the predictions of 
the minimum-energy TIS model. The results of this exercise 
give us some insight into the properties of these halos. In 
Figure 6 we plot the relationship between (t and the mass 
for halos in the VOID and DOUBLE region (the behaviour 
of halos in the SINGLE region is similar to that of those in 
DOUBLE). The difference between halo properties in these 
two regions is striking. In the DOUBLE region there is no 
clear relationship between (^t and mass, but we do not either 
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Figure 8. Plot of ^Qt)- We start plotting from C, = 0.97, where 
the function has an absolute maximuni. The vertical scale is the 
same as that adopted in Fig. 7. 



Figure 9. Measured values of CMo/rtcr^ versus Mass for the 
region of the DOUBLE cluster. Groups marked with a star lie 
within the allowed region where solutions of eq. 13 exist, and the 
corresponding value of C,t can be found. 



find a clustering around the value C,t = 29.4, characterizing 
the minimum- energy TIS solution as suggested by Shapiro 
et al. (1999). On the other hand, there seems to be a rela- 
tionship between (^t and mass for halos in the VOID region, 
although with a rather large dispersion, particularly for ha- 
los having Mo ^ 3 x lO^^/i"^M0. 

A very interesting property of halos in Figure 6 is that 
we do not find halos having (^t > 34.2, the upper limit for 
gravothermal instability for TIS halos, although our upper 
limit for the ^t values extends up to ^t ^ 59.5. We do how- 
ever also find a few halos having ("t < 4.738, the critical value 
below which the total energy E = T + K > {) and the TIS 
solution can not exist (Shapiro et al. 1999). It is important 
to remember that the halos we find in a simulation are not 
spherical and not even "ideal", being a discrete realization 
of some equilibrium state, so the above quoted bounds can 
not be taken literally. 

At first sight, it may seem curious that only a fraction of 
all the halos (46% in the DOUBLE and 37% in the VOID, 
respectively) have values of Gv and Mo for which eq. 13 can 
be solved. The obvious interpretation is that only a frac- 
tion of halos have reached equilibrium, even at the end of 
the simulation; but it also remains the possibility that these 
"out-of-equilibrium" halos relax to an equilibrium state dif- 
ferent from any of the three considered in the present work. 

Finally, in Figs. 10-11 we plot the Gv — Mo relation of 
those halos for which a Ct can be found, for the DOUBLE 
and VOID regions, respectively. Note that for these halos 
the intrinsic dispersion is even smaller than in Figs. 3 and 
4. These halos can be then regarded as very near to a TIS 



equilibrium state. Note also that the coefficient of eq. 1 for 
the TIS solution has been computed for the minimum energy 
TIS solution. In fact, from eq. (98) of Shapiro et al. (1999) 
we see that this coefficient would depend on ^i 



2 
CT„ — 



.2/3 



(37rG)^/^ 0(0 ^2/3 



a(C) 



h;'\i + zcou)m; 



2/3 



(14) 



where: a(C) = 3M(C)/C^p(C)- As we can see from the figures, 
the scatter induced by this dependence is very small and less 
than the intrinsic Poisson scatter. 



4.3 Probability distribution of the spin parameter 

The angular momentum distribution is an interesting stati- 
tistics, because the halo angular momentum originates from 
gravitational interactions between the collapsing region and 
its environment. Following Peebles (1971) and Efstathiou & 
Jones (1979) we will present results for the distribution of 
the spin parameter A defined as: 



A: 



L\E 



1/2 



CMVs 



(15) 



where L and \E\ are the angular momentum and the total 
energy of each halo, respectively. The calculation of |^| is 
not free of ambiguities, because in order to compute the 
potential energy W one should take into account the fact 
that the halo is not isolated, i.e. one should also account for 
the contribution from the environmental gravitational field, 
and this is not currently done by any of the group finders we 
have adopted. For this reason, we show in Figure 12 the spin 



© 2000 RAS, MNRAS 000, 1-13 



10 V. Antonuccio-Delogu et al. 







1 1 1 1 11 1 1 1 1 1 1 1 1 1 , 1 




- 


1000 


- 


DOUBLE - GROUPS 


WITH 


Ct, 






-Shapiro et al. (1999) 




- 






-Bond k Myers (1996) 












: 


SUS 


^^^^.^ 




100 


r 






- 




>^ 


^^ 










I 






10 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 




- 







- 




VOID - GROUPS WITH 


ft , 


1000 






- Shapiro et al. (1999) 


- 










- ooiia (X iviyeib yiuiyo) 




SUS 


^,,r<^- 




8 ^^'-'^ 


_ 




oo JU^'^^""'^ -- 










100 


J^^^^^ .-'■'-" 


- 










,.-■■" oj^^^ 












,--"o"^^/^r^ . " 






.-^^^^ 






I 


' 


10 


- 


- 




- 


- 



1011 5x1011 1012 5x1012 1013 5x1013 

M (M^) 

erouD ^ 0^ 



IQii 5x1011 1012 5x1012 1013 5x1013 

M (M^) 

erouD ^ 0^ 



Figure 10. The dv — Mq relation for groups in the DOUBLE 
region for which the truncation radius can be computed by in- 
verting eq. 13. The scatter is less than in the analogous plot for 
all the groups found by SKID in the DOUBLE region, Fig. 3. 
The error bar in the lower left part of the plot represents the 
uncertainty in the zero point of the best-fit curve induced by the 
scatter in observed values of Qt- Notice that it is less than the 
observed scatter. 



probability distribution P(A) computed only for TIS halos in 
the three regions, i.e. for those halos verifying eq. 13. As we 
have seen in the preceding paragraphs, these halos seem to 
verify the Gv — Mq relationship with a much smaller scatter 
than halos selected by any group finder, so we regard them 
as our fiducial equilibrium halos. Combining eqs. 44 and 45 
from Shapiro et al. (1999) we find that for a TIS halo the 
total energy E is connected to the potential W by: 



2(a- 



-W 



(16) 



The potential W for these halos is then computed exactly, 
i.e. by summing the contribution from each particle in the 
simulation. Note that the parameter a depends on the di- 
mensionless truncation radius Ct which can be evaluated 
only for TIS halos. 

One immediately notices that P(A) seems to depend on the 
environment. It has been shown in recent work that a very 
good fit to P(A) is given by a lognormal distribution (Dal- 
canton et al. 1997; Mo et al. 1998): 



AP(A) 



2llCF\ 



exp 



f ln^(A/(A) l 



d\ 



(17) 



Mo et al. (1998) suggest that Eq. 17 with (A) = 0.05, ax = 
0.05 provides a good fit to the probability distribution of all 
halos, independently of the environment. What we observe 
is, on the other hand, that all the three observed distribu- 



Figure 11. Same as Fig. 10 for the VOID region. Note that the 
scale of the axes is the same as in Figure 5. 



tions seem to be reasonably well fitted by eq. 17, but the 
values of the fitting parameters are certainly very different 
from those mentioned above. Moreover, for those halos se- 
lected with AFOF, the distributions seem to be all consis- 
tent with each other but not well fitted by the lognormal 
model (Figure 13). 

This discrepancy may appear puzzling. However, Figures 12 
and 13 cannot be compared, because of two reasons. First, 
the total energy E entering the definition of A was evaluated 
directly in Figure 12, while it was estimated from eq. 16 in 
Figure 8. The procedure we have adopted to extract fidu- 
cial TIS halos does actually produce a sample which verifies 
the cFv — Mq relationship with a much less statitistical noise 
than the parent sample. But there is an even stronger reason 
which makes the comparison doubtful: the total number of 
halos extracted using AFOF is much less than that obtained 
using SKID. Moreover, the extent over A of the spin prob- 
ability distribution is smaller than for SKID as is evident 
from a comparison of the figures. If we take these differences 
into account, we do not see any significant difference among 
the distributions in the VOID region. For all these reasons, 
we can conclude that there is a dependence of the spin prob- 
ability distribution P(A) on the environment only for TIS 
halos. It would be interesting to speculate about the physi- 
cal mechanisms producing this dependence, and we hope to 
be able to address this question in further work. 

Before closing this section, we would like to remind the 
reader that recent theoretical calculations predict a rather 
large distribution in the average values and shape of P(A) , 
with a rather marked dependence on the peak's overdensity 
(Catelan & Theuns 1996) or on the details of the merging 
histories (Nagashima & Gouda 1998; Vitvitska et al. 2001). 
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Table 3. Fits of P(A) with lognormal distribution, for TIS halos. 
Region Median (A) ax 

VOID 0.018 0.018 0.5 

DOUBLE 0.018 0.03 0.9 

SINGLE 0.051 0.07 1.4 

All (SKID) 0.06 0.07 0.65 
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Figure 12. Probability distribution of the spin parameter for the 
three regions, for TIS halos. The dashed curves are the best-fit ap- 
proximations obtained using the lognormal distribution adopted 
by Mo et al. (1998) (eq. 15). Histograms and fitting curves are 
normalized to unity. 



A direct comparison of our results with the conditional prob- 
ability distribution P(A|z/) of Catelan & Theuns (1996) is 
made difficult by the fact that the relationship between the 
linear overdensity v and (for instance) the final mass of the 
halo turns out to be quite noisy (Sugerman et al. 2000, Fig. 
10), so it is not possible to "label" unambigously each halo 
with its initial overdensity. However, one could hope to fur- 
ther increase the number of halos by further diminishing the 
softening length, and we hope to get a better statistics from 
future simulations which would help us to address also the 
latter points. 



4.4 Density profiles of massive halos 

As we already mentioned in the Introduction, even the most 
massive halos we find in this simulation using SKID do not 
contain enough particles to allow a reliable determination 
of the density profile. This is clearly visible from Figure 13, 
where we plot the profiles of the four most massive halos 
extracted from the DOUBLE cluster region. None of these 
halos lies in the integrability strip of Figure 11, so the best 
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Figure 13. Same as Figure 14 but for halos selected using AFOF. 
The dashed curve is a lognormal distribution with }A(= 0.5, ay = 
0.05. 



fit TIS profiles displayed as continous curves have been ob- 
tained by least square fittings, where we have varied po and 
ro. 

The distinguishing feature of the TIS density profile, 
when compared with the universal density profile of Navarro 
et al. (1996, 1997) is the presence of a central core. Although 
the minimum energy TIS profiles fit reasonably well the cen- 
tral regions, they fall out too gently at distances larger than 
a few times the core radius, and in no case we can find a 
reasonably good overall agreement. It would be hazardous 
to draw any conclusion from this comparison, in view of the 
above mentioned poor resolution. However, a reasonable ex- 
planation for the sharp decline of the density profiles is tidal 
stripping, which should be effective at a few times the core 
radius. 



5 CONCLUSION 

The properties of galaxy-sized halos we have considered in 
this paper seem to be very constraining for halo collapse 
and equilibrium models. However, none of the equilibrium 
models considered (neither the minimum energy TIS model) 
seems to be able to give a comprehensive description of our 
findings. We would like to summarize now our findings and 
to point to some controversial issues that they pose. 
First of all, the av — Mq statistics seems to be a sensitive 
tool to discriminate among different halo equilibrium mod- 
els. This statistics is easy to evaluate, because it relies on 
global quantities, and it can then be applied to samples of 
halos. In this context, it is more difficult to discriminate 
models using statistics like the density profile, which would 
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Figure 14. Density profiles for the four most massive halos 
within the DOUBLE region. The number of particles contained in 
each halo is shown. Best fit solutions using the minimum energy 
TIS solution by Shapiro et al. (1999) are shown. 



require a considerably larger mass range in order to give re- 
liable results (see for instance Jing & Suto 2000). 
Models for the Gv — Mq statistics based on the NFW density 
profile seem to be only marginally consistent with simulation 
data. The role of the anisotropy parameter in this context 
does not seem to be crucial: it is the slope of the radius-mass 
relationship for these halos which seems to mostly affect the 
normalisation of the av — Mq . 

As we have seen, the TIS model seems to offer a very good 
quantitative framework to explain the ay — Mq statistics, 
even in the VOID region where the slope of the relation- 
ship is very different from that predicted by the minimum 
energy TIS model of Shapiro et al. The fact that a model 
based on the hypothesis that halos have a finite extent pro- 
vides a good description should not come as a surprise. Halos 
forming in clusters experience a complex tidal field originat- 
ing from neighbouring halos and from the large scale web 
in which they are embedded. The tidal radii of the environ- 
ments within which they lie, although often larger than the 
mean distance, could limit the extent of halos. A theoret- 
ical treatment of the growth of the angular momentum is 
complicated by the fact that the distribution of the torques 
induced by nearby halos depends on clustering (Antonuccio- 
Delogu & Atrio-Barandela 1992). However, we believe that 
it would be difficult to think that the truncation is a nu- 
merical artifact due to the finite mass resolution: were this 
the case, we should expect the same relationship between 
truncation radius rt and mass in all the three regions, but 
this is clearly not the case. 
We have already noted the fact that the truncation radii 



we find are always less than the critical value for the on- 
set of gravot hernial instability, Ccrit = 34.2. This leads us 
to think that this instability is at work in our simulations, 
but in order to investigate this issue one would need simula- 
tions with a dynamical range at least 3 orders of magnitudes 
larger than those used in this simulation. 
Concerning the dependence of the distribution of spin pa- 
rameters on the environment, we find that halos selected 
using A EOF do not show any dependence on the environ- 
ment (the same holds for halos selected using SKID), but 
if we select subsamples of fiducial TIS halos, we do find a 
dependence of the properties of P(A) on the environment. In 
particular, this fact seems to be at odds with the recent in- 
vestigation by Syer et al. (1999), who find that the observed 
distribution of the spin parameter for a large homogeneous 
sample of spirals is well described by a lognormal distribu- 
tion with (A) ~ 0.05 and a variance ax ~ 0.36. This result is 
in contrast also with other work mentioned above (Warren 
et al. 1992; Eisenstein & Loeb 1995). If confirmed by further 
investigations, this discrepancy could suggest that there is 
probably some systematic trend in the way the angular mo- 
mentum of the luminous discs is connected with that of the 
halo, which is not accounted for by the models of Syer et al. 
(1999). 

Last but not least, it is important to stress that Lemson & 
Kauffmann (1999) conclude that "Only the mass distribu- 
tion varies as a function of environment. This variation is 
well described by a simple analytic formula based on the 
conditional Press-Schechter theory. We find no significant 
dependence of any other halo property on environment...". 
In comparing their results with ours, we must keep in mind 
that we have followed a very different procedure from theirs, 
because we have prepared a simulation using constrained ini- 
tial conditions with the purpose of obtaining a final configu- 
ration containing certain features (i.e. a double cluster and 
a void). Although our simulation box is not an "average" 
region of the Universe, it is certainly a representative one. 
We stress again the fact that all the halos from underdense 
regions in our simulation come from a void, and not from 
the outer parts of clusters. Lemson & Kauffmann, on the 
other hand, seem to take their halos from all the volume 
and group them according to the overdensity of their parent 
regions. We think then that a direct comparison between 
the results of these two different investigations would be 
misleading, given the complementarity of our approaches. 



6 APPENDIX 

We give a full derivation of equation 13. The starting point 
is the virial theorem for systems with boundary pressure 
terms, as given in eq. 41 from Shapiro et al. (1999): 



= 2K + W + Sp 



(18) 



In the above equation the kinetic energy K can be rewritten 
in terms of the 1-D velocity dispersion: 



Mojvf 3 2 

= 2 "" 2 °^ 

The potential energy term W: 



(19) 
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W = AnG I pM{r)rdr = AtvG / pM{r)rdr (20) 

can be rewritten in terms of global quantities and of the 
dimensionless radius (t : 

GMi Ct*(0) 



W = 

where we have defined: 

*(C0 = / ' dCCp(C)M(C) 



(21) 



Finally, Sp is a surface term which arises from the constraint 
that the system has a finite radius, and is given by (Shapiro 
et al., eq. 43): 



Sp = -AnroptCt 



(22) 



We are adopting here the same notation as Shapiro et al., 
so that ro and pt are the core radius and an external "pres- 
sure" term, respectively. Using eqs. 34 and 38 from Shapiro 
et al., the latter equation can be rewritten in terms of the 
dimensionless integrated mass and density: 

^P = -^i^TTT^^^*)'^' (23) 

Substituting eqs. 19, 21, 23 into eq. 18 we get: 

^A^^2 GMJ Ct^(Ct) Mo 2 f^,. 

from which we get the desired equation. 
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